#include <iostream>
#include <vector>
#include <list>
#include <functional>
#include <numeric>
#include <TClonesArray.h>
#include <TRefArray.h>
#include <TMath.h>
#include <TF1.h>

#include "TROOT.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TChain.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TString.h"
#include <TLegend.h>
#include "TLine.h"
#include "TH1F.h"
#include "TH2F.h"
#include "bunchtime.cxx"
#include "fluxweight.cxx"

using namespace std;
#if !defined(__CINT__) || defined(__MAKECINT__)
#include "AnalysisTools/TCCreconAnalysisUtils.hxx"
#include "AnalysisTools/TGlobalAnalysisUtils.hxx"
#include "AnalysisTools/TNEUTAnalysisUtils.hxx"
#include "AnalysisTools/TTrackerAnalysisUtils.hxx"
#include "AnalysisTools/TCCtruthAnalysisUtils.hxx"
#include "AnalysisTools/TCheatingAnalysisUtils.hxx"
#include "AnalysisTools/TCheatingAnalysisUtils.hxx"
#include "AnalysisTools/TGlobalAnalysisUtils.hxx"
#include "AnalysisTools/TTruthUtils.hxx"
#endif

using std::cout;
using std::endl;

//This is the start of the analysis program.  Leave everything above
//here alone.  So the name after void is the name that should be the
//same as the macro name.This is the name of the thing you run in
//root after running .x oaAnalysisInitFile.C and then Macro.C
//Filelist is the name of the .txt or .dat file that contains the
//list of MC or data files you want to analyse.
void mymacro12_edit(char *filelist){
    
    //*****************************
    // Define FGD fiducial volumes (all dimensions in mm) we'll use these numbers in cuts later:
    //*****************************
    const float FGD1xmin = -874.51;
    const float FGD1ymin = -819.51;
    const float FGD1zmin = 136.875;
    const float FGD1xmax = 874.51;
    const float FGD1ymax = 929.51;
    const float FGD1zmax = 446.955;
    float FGD1xmid, FGD1ymid;
    FGD1xmid = (FGD1xmin+FGD1xmax)/2.0;
    FGD1ymid = (FGD1ymin+FGD1ymax)/2.0;
    
    //*************************************
    // Define constants for Enu calculation at end
    //*************************************
    //energy reconstruction and direction of lepton wrt beam for 1 ring events
    const double V  = 27.0;            // MeV,  binding potential ********** NEEDS CHANGING TO 25 WHEN CHECKS DONE ***********
    const double mn = 939.56536;       // neutron mass
    const double mp = 938.27203;       // proton mass
    const double ml = 105.65837;       // muon mass
    // const double ml = 0.510999;        // electron mass
    // const double nudir[3] = {0.5486540, -0.835246, 0.0366437}; //beam direction - SK1
    // const double nudir[3] = {0.669764, -0.742179, 0.024223}; //beam direction - SK Tom
    const double nudir[3] = {0.0, 0.0, 1.0}; //beam direction - nd280, assuming along z axis
    const double mumass = 105.65837;   // muon mass
    const double pimass = 139.57018;   // pi ch mass
    
    //*******************
    // Define Histograms
    //*******************
    
    const int nbins = 15;
    double binEdges[nbins+1]={0,400,800,1200,1600,2000,2800,3600,4400,5200,6000,6800,7600,8400,9200,10000};
    
    TH1F *hEnu    = new TH1F("hEnu", "Energy est CCQE #nu_{#mu}", 1000, 0, 10000);
    TH1F *hEest_all = new TH1F("hEest_all", "Energy est all #nu_{#mu}", 1000, 0, 10000);  //10000
    TH1F *hpbunch = new TH1F("hpbunch", "Energy est nonCCQE #nu_{#mu}", 1000, 0, 10000); //10000
    
    TH1F *hEestQ1 = new TH1F("hEestQ1", "Energy est all #nu_{#mu}", 1000, 0, 10000); //MC version of data histos
    TH1F *hEestQ3 = new TH1F("hEestQ3", "Energy est all #nu_{#mu}", 1000, 0, 10000); //MC version of data histos
    TH1F *hEestQ1_b = new TH1F("hEestQ1_b", "Energy est all #nu_{#mu}", nbins, binEdges); // same but binned
    TH1F *hEestQ3_b = new TH1F("hEestQ3_b", "Energy est all #nu_{#mu}", nbins, binEdges); // same but binned
    
    TH2F *hEnuVsTrue = new TH2F("hEnuVsTrue", "CCQE Enu test", 1000, 0, 10000, 1000, 0, 10000);
    TH2F *hpbunchVsTrue = new TH2F("hpbunchVsTrue", "nonCCQE - vecsum test", 1000, 0, 10000, 1000, 0, 10000);
    TH2F *hpbunchcompVsTrue = new TH2F("hpbunchcompVsTrue", "nonCCQE - vecsum + masses test", 1000, 0, 10000, 1000, 0, 10000);
    TH2F *hpsumcompVsTrue = new TH2F("hpsumcompVsTrue", "nonCCQE - scalsum + masses test", 1000, 0, 10000, 1000, 0, 10000);
    TH2F *hEtracksumVsTrue = new TH2F("hEtracksumVsTrue", "nonCCQE - sum E of tracks", 1000, 0, 10000, 1000, 0, 10000);
    
    TH1F *hreconspecPiQ1 = new TH1F("hreconspecPiQ1", "Recon energy per parent per quad", 1000, 0, 10000);
    TH1F *hreconspecPiQ3 = new TH1F("hreconspecPiQ3", "Recon energy per parent per quad", 1000, 0, 10000);
    TH1F *hreconspecKQ1 = new TH1F("hreconspecKQ1", "Recon energy per parent per quad", 1000, 0, 10000);
    TH1F *hreconspecKQ3 = new TH1F("hreconspecKQ3", "Recon energy per parent per quad", 1000, 0, 10000);
    TH1F *hreconspecPiQ1_b = new TH1F("hreconspecPiQ1_b", "Recon energy per parent per quad", nbins, binEdges);
    TH1F *hreconspecPiQ3_b = new TH1F("hreconspecPiQ3_b", "Recon energy per parent per quad", nbins, binEdges);
    TH1F *hreconspecKQ1_b = new TH1F("hreconspecKQ1_b", "Recon energy per parent per quad", nbins, binEdges);
    TH1F *hreconspecKQ3_b = new TH1F("hreconspecKQ3_b", "Recon energy per parent per quad", nbins, binEdges);
    TH1F *hreconspecMuQ1 = new TH1F("hreconspecMuQ1", "Recon energy per parent per quad", 1000, 0, 10000);// just checking contribution
    TH1F *hreconspecMuQ3 = new TH1F("hreconspecMuQ3", "Recon energy per parent per quad", 1000, 0, 10000);// just checking contribution
    
    
    hEnu ->SetLineColor(kRed+1);
    hEnu ->GetXaxis()->SetTitle("E est (MeV) for CCQE #nu");
    hEnu ->GetYaxis()->SetTitle("#");
    hpbunch ->SetLineColor(kBlue+1);
    hpbunch ->GetXaxis()->SetTitle("E est (MeV)for non-CCQE #nu");
    hpbunch ->GetYaxis()->SetTitle("#");
    
    hEestQ1 ->SetLineColor(kBlue+1);
    hEestQ1 ->GetXaxis()->SetTitle("E est (MeV)for #nu_{#mu}");
    hEestQ1 ->GetYaxis()->SetTitle("#");
    hEestQ3 ->SetLineColor(kRed+1);
    hEestQ3 ->GetXaxis()->SetTitle("E est (MeV)for #nu_{#mu}");
    hEestQ3 ->GetYaxis()->SetTitle("#");
    hEestQ1_b ->SetLineColor(kBlue+1);
    hEestQ1_b ->GetXaxis()->SetTitle("E est (MeV)for #nu_{#mu}");
    hEestQ1_b ->GetYaxis()->SetTitle("#");
    hEestQ3_b ->SetLineColor(kRed+1);
    hEestQ3_b ->GetXaxis()->SetTitle("E est (MeV)for #nu_{#mu}");
    hEestQ3_b ->GetYaxis()->SetTitle("#");
    
    hEnuVsTrue ->GetXaxis()->SetTitle("True E (MeV)");
    hEnuVsTrue ->GetYaxis()->SetTitle("Enu");
    hpbunchVsTrue ->GetXaxis()->SetTitle("True E (MeV)");
    hpbunchVsTrue ->GetYaxis()->SetTitle("pbunch");
    hpbunchcompVsTrue ->GetXaxis()->SetTitle("True E (MeV)");
    hpbunchcompVsTrue ->GetYaxis()->SetTitle("vecmom+masses");
    hpsumcompVsTrue ->GetXaxis()->SetTitle("True E (MeV)");
    hpsumcompVsTrue ->GetYaxis()->SetTitle("scalmom+masses");
    hEtracksumVsTrue ->GetXaxis()->SetTitle("True E (MeV)");
    hEtracksumVsTrue ->GetYaxis()->SetTitle("sum E tracks");
    
    hreconspecPiQ1 ->SetLineColor(kBlue+1);
    hreconspecPiQ1 ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecPiQ1 ->GetYaxis()->SetTitle("#");
    hreconspecPiQ3 ->SetLineColor(kRed+1);
    hreconspecPiQ3 ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecPiQ3 ->GetYaxis()->SetTitle("#");
    hreconspecKQ1 ->SetLineColor(kBlue+1);
    hreconspecKQ1 ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecKQ1 ->GetYaxis()->SetTitle("#");
    hreconspecKQ3 ->SetLineColor(kRed+1);
    hreconspecKQ3 ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecKQ3 ->GetYaxis()->SetTitle("#");
    hreconspecPiQ1_b ->SetLineColor(kBlue+1);
    hreconspecPiQ1_b ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecPiQ1_b ->GetYaxis()->SetTitle("#");
    hreconspecPiQ3_b ->SetLineColor(kRed+1);
    hreconspecPiQ3_b ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecPiQ3_b ->GetYaxis()->SetTitle("#");
    hreconspecKQ1_b ->SetLineColor(kBlue+1);
    hreconspecKQ1_b ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecKQ1_b ->GetYaxis()->SetTitle("#");
    hreconspecKQ3_b ->SetLineColor(kRed+1);
    hreconspecKQ3_b ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecKQ3_b ->GetYaxis()->SetTitle("#");
    hreconspecMuQ1 ->SetLineColor(kBlue+1);
    hreconspecMuQ1 ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecMuQ1 ->GetYaxis()->SetTitle("#");
    hreconspecMuQ3 ->SetLineColor(kRed+1);
    hreconspecMuQ3 ->GetXaxis()->SetTitle("E recon (MeV)for #nu_{#mu}");
    hreconspecMuQ3 ->GetYaxis()->SetTitle("#");
    
    
    TLegend* legemuQ1Q3;
    legemuQ1Q3 = new TLegend(0.6,0.7,0.89,0.89);  // fractions of pad dimensions
    legemuQ1Q3->AddEntry(hEestQ1,"all #nu_{#mu} in Q1","l");  // "l" means line "f" for box
    legemuQ1Q3->AddEntry(hEestQ3,"all #nu_{#mu} in Q3","l");
    legemuQ1Q3->SetHeader("Spectra of all #nu_{#mu} across FGD1");
    
    TLegend* legePiKQ1Q3;
    legePiKQ1Q3 = new TLegend(0.6,0.7,0.89,0.89);  // fractions of pad dimensions
    legePiKQ1Q3->AddEntry(hreconspecPiQ1,"#nu_{#mu} from #pi^{#pm} in Q1","l");  // "l" means line "f" for box
    legePiKQ1Q3->AddEntry(hreconspecPiQ3,"#nu_{#mu} from #pi^{#pm} in Q3","l");
    legePiKQ1Q3->AddEntry(hreconspecKQ1,"#nu_{#mu} from K(any) in Q1","l");
    legePiKQ1Q3->AddEntry(hreconspecKQ3,"#nu_{#mu} from K(any) in Q3","l");
    legePiKQ1Q3->SetHeader("Spectra of #nu_{#mu} across FGD1 by parent");
    TLegend* legePiKQ1Q3_b;
    legePiKQ1Q3_b = new TLegend(0.6,0.7,0.89,0.89);  // fractions of pad dimensions
    legePiKQ1Q3_b->AddEntry(hreconspecPiQ1_b,"#nu_{#mu} from #pi^{#pm} in Q1","l");  // "l" means line "f" for box
    legePiKQ1Q3_b->AddEntry(hreconspecPiQ3_b,"#nu_{#mu} from #pi^{#pm} in Q3","l");
    legePiKQ1Q3_b->AddEntry(hreconspecKQ1_b,"#nu_{#mu} from K(any) in Q1","l");
    legePiKQ1Q3_b->AddEntry(hreconspecKQ3_b,"#nu_{#mu} from K(any) in Q3","l");
    legePiKQ1Q3_b->SetHeader("Spectra of #nu_{#mu} across FGD1 by parent");
    
    //*****************
    // Define canvases
    //*****************
    
    TCanvas *c1  = new TCanvas("c1","c1",500,250);
    TCanvas *c2  = new TCanvas("c2","c2",500,250);
    TCanvas *c3  = new TCanvas("c3","c3",500,250);
    TCanvas *c4  = new TCanvas("c4","c4",500,250);
    TCanvas *c5  = new TCanvas("c5","c5",500,250);
    TCanvas *c6  = new TCanvas("c6","c6",500,250);
    TCanvas *c7  = new TCanvas("c7","c7",500,250);
    TCanvas *c8  = new TCanvas("c8","c8",500,250);
    TCanvas *c9  = new TCanvas("c9","c9",500,250);
    TCanvas *c10  = new TCanvas("c10","c10",500,250);
    TCanvas *c11  = new TCanvas("c11","c11",500,250);
    TCanvas *c12  = new TCanvas("c12","c12",500,250);

    
    //**********************
    //Prepare to read trees.
    //**********************
    //
    // Make a chain for the Tracker tree stored in the ReconstructionDir. This needs
    // to be done for each analysis tree you want to get information from.
    //
    //Version for non-skimmed files
    //TChain reconTracker         = TChain("ReconDir/Tracker");
    //TChain mcTruthRooTrackerVtx = TChain("TruthDir/NRooTrackerVtx");    // GDB [GN]
    //TChain mcTruthVertices      = TChain("TruthDir/Vertices");          //HMOK
    
    //Using bits of tracker and global recon - global generally better but tracker better for 1 bit...
    TChain reconTracker     =TChain("ReconDir/Tracker");
    TChain reconGlobal      = TChain("ReconDir/Global");             //HMOK
    TChain BasicHeader      = TChain("HeaderDir/BasicHeader");
    TChain BasicDataQuality = TChain("HeaderDir/BasicDataQuality");
    TChain BeamSummaryData  = TChain("HeaderDir/BeamSummaryData");
    
    TChain mcTruthRooTrackerVtx = TChain("TruthDir/NRooTrackerVtx");    // GDB [GN]
    TChain mcTruthVertices      = TChain("TruthDir/Vertices");          //HMOK
    TChain mcTruthTrajects      = TChain("TruthDir/Trajectories");      //HMOK
    
    //Sets up the pointer to the bunch calculation class/function thing
    //Calculates bunch based on time cos bunch number not given
    bunchtime   buntm;
    fluxweight  fluxw;
    
    //Here is where we read in the list of files from the filelist then add branches you need info from
    int k = 0;  //counter for number of files.
    char filename[128]; //This sets up a variable to contain the text in filelist [128] is related to the length of text
    string prevfilename;
    ifstream in;
    in.open(filelist);
    while (in){
        in >> filename;
        if(k == 0){
            reconTracker.AddFile(filename);
            //Other branches you wish to chain would go in here.
            reconGlobal.AddFile(filename);        // HMOK Added for recon versions
            BasicHeader.AddFile(filename);
            BasicDataQuality.AddFile(filename);
            BeamSummaryData.AddFile(filename);
            mcTruthRooTrackerVtx.AddFile(filename);   // GDB
            mcTruthVertices.AddFile(filename);        // HMOK
            mcTruthTrajects.AddFile(filename);
            prevfilename = filename;
        }
        
        //This next piece of logic makes sure we don't read the last file in twice.
        if(k != 0 && prevfilename != filename){
            reconTracker.AddFile(filename);
            //Other branches you wish to chain would go in here.
            BasicHeader.AddFile(filename);
            BasicDataQuality.AddFile(filename);
            BeamSummaryData.AddFile(filename);
            reconGlobal.AddFile(filename);        // HMOK
            mcTruthRooTrackerVtx.AddFile(filename);   // GDB
            mcTruthVertices.AddFile(filename);        // HMOK
            mcTruthTrajects.AddFile(filename);

        }
        prevfilename = filename;
        ++k;
    }
    in.close();
    
    
    //This sets up a TClonesArray to allow us to access the information from the
    //the branchs we chained above.
    
    Int_t EventID;
    Int_t RunID;
    UInt_t EventTime;
    BasicHeader.SetBranchAddress("EventTime", &EventTime);
    BasicHeader.SetBranchAddress("EventID", &EventID);
    BasicHeader.SetBranchAddress("RunID", &RunID);
    
    Int_t DQEventID;
    Int_t DQRunID;
    Int_t ND280OffFlag;
    BasicDataQuality.SetBranchAddress("EventID", &DQEventID);
    BasicDataQuality.SetBranchAddress("RunID", &DQRunID);
    BasicDataQuality.SetBranchAddress("ND280OffFlag", &ND280OffFlag);
    
    bool health = true; //no idea what this is for...
    Int_t ND280Spill;
    Int_t BRunID;
    Int_t BEventID;
    Int_t BeamSummaryDataStatus;
    
    BeamSummaryData.SetBranchAddress("RunID",&BRunID);
    BeamSummaryData.SetBranchAddress("EventID",&BEventID);
    BeamSummaryData.SetBranchAddress("BeamSummaryDataStatus",&BeamSummaryDataStatus);
    BeamSummaryData.SetBranchAddress("ND280Spill", &ND280Spill);
    
    TClonesArray *TClonesBSD = new TClonesArray("ND::TBeamSummaryDataModule::TBeamSummaryData", 20);
    BeamSummaryData.SetBranchAddress("BeamSummaryData",&TClonesBSD);
    ND::TBeamSummaryDataModule::TBeamSummaryData *bsd=(ND::TBeamSummaryDataModule::TBeamSummaryData*)(*TClonesBSD)[0];
    
    //NTracks matches the name in the oaAnalysis .root files and is the number of tracks in the tracker per MC entry
    Int_t NPIDs, NVertices;
    TClonesArray* PIDs = new TClonesArray("ND::TGlobalReconModule::TGlobalPID", 200);
    TClonesArray* GVertices = new TClonesArray("ND::TGlobalReconModule::TGlobalVertex", 200);
    reconGlobal.SetBranchAddress("NVertices",&NVertices);
    reconGlobal.SetBranchAddress("Vertices",&GVertices);
    reconGlobal.SetBranchAddress("NPIDs",&NPIDs);
    reconGlobal.SetBranchAddress("PIDs",&PIDs);
    
    Int_t NTracks;
    TClonesArray *Tracks = new TClonesArray("ND::TTrackerReconModule::TTrackerResult", 20);
    reconTracker.SetBranchAddress("NTracks", &NTracks);
    reconTracker.SetBranchAddress("Tracks", &Tracks);
    
    Int_t NVtx;  // GDB
    TClonesArray *Vtx = new TClonesArray("ND::NRooTrackerVtx", 100);  // GDB [GN]
    if (Vtx == 0) cout << "p_TruthRoos TClonesArray pointer is null\n";
    mcTruthRooTrackerVtx.SetBranchAddress("NVtx",&NVtx);   // GDB
    mcTruthRooTrackerVtx.SetBranchAddress("Vtx",&Vtx);  // GDB
    
    Int_t NTraj;
    TClonesArray *Trajectories = new TClonesArray("ND::TTruthTrajectoriesModule::TTruthTrajectory", 20);
    mcTruthTrajects.SetBranchAddress("NTraj", &NTraj);
    mcTruthTrajects.SetBranchAddress("Trajectories", &Trajectories);
    
    Int_t NTVtx;
    TClonesArray *Vertices = new TClonesArray("ND::TTruthVerticesModule::TTruthVertex", 20);
    mcTruthVertices.SetBranchAddress("NVtx", &NTVtx);
    mcTruthVertices.SetBranchAddress("Vertices", &Vertices);
    
    //Define stuff used later in loops that I might want to plot @ end
    int tpcno;
    float  tpcmom, tpcchr;
    
    //Gets the number of entries in the chained file (sum of entries for all files used in filelist)
    int nentries = reconTracker.GetEntries();
    std::cout << "Number entries: " << nentries << std::endl;
    //###### Here! ######
    //nentries = 5000; //use this to set number of entries if want limited test run
    int nTracker = 0;
    int ntpcpassall = 0, npasspulle = 0, npasstpcmom = 0, npasstracks = 0;
    int reconmuon = 0, truemuon = 0;
    //counters
    int Nnomu = 0, NCCQEmu = 0, NnCCQEmu = 0;
    
    
    BasicHeader.GetEntry(0);
    BasicDataQuality.GetEntry(0);
    BeamSummaryData.GetEntry(0);
    mcTruthRooTrackerVtx.GetEntry(0);  // GDB
    mcTruthVertices.GetEntry(0);  // GDB
    mcTruthTrajects.GetEntry(0); //HMOK
    int entries = BasicHeader.GetEntries();
    
    int EventsSeen = 0;
    int EventsPassedPODVeto = 0;
    int EventsPassedTPC1Veto = 0;
    int EventsPassedBunchTiming = 0;
    int EventsPassedTPCInfo = 0;
    int EventsPassedVertexCut = 0;
    int EventsPassedForwardCut = 0;
    int EventsPassedTPC2Comp = 0;
    int EventsPassedMuonCandCut = 0;
    int EventsPassedMuonCut = 0;
    int EventsPassedCCQE = 0;
    int EventsPassednonCCQE = 0;
    int EmptyBunch = 0;
    
    int bob = 1304156;
    
    const float bunch1 = 2750.2;
    const float bunch2 = 3332.0;
    const float bunch3 = 3914.7;
    const float bunch4 = 4497.0;
    const float bunch5 = 5078.4;
    const float bunch6 = 5659.7;
    const float bunch7 = 6243.4;
    const float bunch8 = 6824.2;

    
    
    
    //Load in the histograms for the flux reweighting
    TFile *flux = new TFile("/home/okeeffe/tuned11bv3.2/run2/nd5_tuned11bv3.2_11anom_run2_fine.root");//beamb
    //TFile *flux = new TFile("/home/okeeffe/tuned11bv3.2/run3c/nd5_tuned11bv3.2_11anom_run3c_fine.root");//beamc
    TH1F *numu  = (TH1F*)flux->Get("enu_nd5_tuned11b_numu_ratio");
    TH1F *numub = (TH1F*)flux->Get("enu_nd5_tuned11b_numub_ratio");
    TH1F *nue   = (TH1F*)flux->Get("enu_nd5_tuned11b_nue_ratio");
    TH1F *nueb  = (TH1F*)flux->Get("enu_nd5_tuned11b_nueb_ratio");

    
    Double_t POTs(0);
    Int_t eveTot(0), eveBeam(0), eveDQ(0);
    //Start running over all entries in the chained file
    for(int ient = 0; ient < nentries; ient++){                   //  *IN ENTRIES LOOP*
        //Add event ticker
        if(ient%10000 == 0) cout << "Processing event number " << ient << endl;
        
        //Access the information in the branchs of interest for entry ient
        reconTracker.GetEntry(ient);
        reconGlobal.GetEntry(ient);
        BasicHeader.GetEntry(ient);
        BasicDataQuality.GetEntry(ient);
        BeamSummaryData.GetEntry(ient);
        mcTruthRooTrackerVtx.GetEntry(ient);  // GDB
        mcTruthVertices.GetEntry(ient);  // GDB
        mcTruthTrajects.GetEntry(ient); //HMOK
//        if (RunID != BRunID || EventID != BEventID || RunID != DQRunID || EventID != DQEventID ) {std::cout << "Runs and Event # mismatch." << std::endl;}
//        //***ASSESS WHETHER GOOD ND280 DATA
//        // Only data quality assessed data (ND280OffFlag = 0 ), Beam Data OK (BeamSummaryDataStatus == 1), and good spill ((Int_t (bsd->GoodSpillFlag)) == 1)
//        if ( ND280OffFlag == 0 && BeamSummaryDataStatus == 1 && (Int_t (bsd->GoodSpillFlag)) == 1 )
//            POTs += bsd->CT5ProtonsPerSpill;
//        if ( ND280OffFlag == 0 && BeamSummaryDataStatus == 1 && (Int_t (bsd->GoodSpillFlag)) == 1 )
//            eveTot += 1;
//        if ( ND280OffFlag == 0 && BeamSummaryDataStatus == 1 && (Int_t (bsd->GoodSpillFlag)) == 1 )
//            eveBeam += 1;
//        if ( ND280OffFlag == 0 && BeamSummaryDataStatus == 1 && (Int_t (bsd->GoodSpillFlag)) == 1 )
//            eveDQ += 1;
//        bool goodevent = ND280OffFlag == 0 && BeamSummaryDataStatus == 1 && (Int_t (bsd->GoodSpillFlag)) == 1;
//        cout << "Checking if good event" << endl;
//        if(!goodevent){
//            cout << "Not good event, skipping" << endl;
//            continue;}
//        else {cout << "Good event, continuing..." << endl;}

               
        //Loop through the PIDs to see how many P0D, P0DECal and TPC1 objects there are in the bunch.
        //We want to veto events where any of these are non-zero.
        int pod[8]  = {0,0,0,0,0,0,0,0};
        int tpc1a[8] = {0,0,0,0,0,0,0,0};
        int pode[8] = {0,0,0,0,0,0,0,0};
        
        for(int bun = 0; bun < 8; ++bun){
            EventsSeen++;
            int peb = 0, tpc1bb = 0, tpc2bb = 0, tpc3bb = 0, pdb = 0;
            for(int nt = 0; nt < NTracks; ++nt){
                //Loop over all tracks. Using trTrack to access Tracks[nt] - each track object?
                ND::TTrackerReconModule::TTrackerResult *trTrack = (ND::TTrackerReconModule::TTrackerResult*)(*Tracks)[nt];
                
                if(!((fabs(trTrack->Position.T() - bunch1) < 60.) || (fabs(trTrack->Position.T() - bunch2) < 60.) || (fabs(trTrack->Position.T() - bunch3) < 60.) || (fabs(trTrack->Position.T() - bunch4) < 60.) || (fabs(trTrack->Position.T() - bunch5) < 60.) || (fabs(trTrack->Position.T() - bunch6) < 60.) || (fabs(trTrack->Position.T() - bunch7) < 60.) || (fabs(trTrack->Position.T() - bunch8) < 60.))){continue;}
                
                
                // Feeding in track time to bunch calculator bunchtime.cxx
                // If bunch no. calculated from track position in time is not same as current bunch continue
                if(bun != buntm.bunchno(trTrack->Position.T())){continue;}
                // Create TClones array TPC of TPC objects for each track?
                TClonesArray * TPC = trTrack->TPC;
                // Loop over number of TPC objects per track
                for(int ntpc = 0; ntpc < trTrack->NTPCs; ntpc++){
                    // tpctrack accesses each TPC object and checks which detector it originated it
                    ND::TTrackerReconModule::TTPCTrack *tpctrack = (ND::TTrackerReconModule::TTPCTrack*)(*TPC)[ntpc];
                    if(tpctrack->Detector == 1){++tpc1bb;}
                    if(tpctrack->Detector == 2){++tpc2bb;}
                    if(tpctrack->Detector == 3){++tpc3bb;}
                    if(EventID==bob){std::cout << "Momentum: " << trTrack->Momentum << " Bunch: " << buntm.bunchno(trTrack->Position.T()) << " Det: " << tpctrack->Detector << std::endl;}
                }
            }
            tpc1a[bun] = tpc1bb; // Stores number of TPC1 objects per bunch
            
            // Loop over objects to describe a reconstructed PID
            for(int ngt = 0; ngt < NPIDs; ++ngt){
                // Gets each PID object
                ND::TGlobalReconModule::TGlobalPID *gbTrack = (ND::TGlobalReconModule::TGlobalPID*)(*PIDs)[ngt];
                // Gets front position of track and checks if timing correct
                if(bun != buntm.bunchno(gbTrack->FrontPosition.T())){
                    continue;}
                // Checks number of POD objects for this PID object
                pdb = pdb+gbTrack->NP0Ds;
                //if(EventID==bob){std::cout << "Momentum: " << gbTrack->FrontMomentum << " Bunch: " << bun << " nPOD: " << pdb << std::endl;}
            }
            pod[bun] = pdb; // Stores number of POD objects per bunch
            
        }
        //***** RECON SECTION *****
        //Now loop over the number of tracks in the event to get information about them
        // and to see if they satisfy the cuts/selection criteria
        
        //Need FGDtiming bin info for Michel CCQE cut, so do second quick loop over everything in entry
        //  for(){
        //    //get bunch peak positions, define boundaries using +-140ns
        //  }
        
        
        for(int bunch = 0; bunch < 8; ++bunch){                       //  **IN BUNCH LOOP**
            if(EventID==bob){std::cout << "Marker 0" << endl;}
            //std::cout << "Event: " << EventID << std::endl;
            
            //Define quantities that will use at end of PID loop
            float pmax_passmu = 0.0, pmax_notmu = 0.0;
            float xpos_passmu, xpos_notmu;
            float ypos_passmu, ypos_notmu;
            float zpos_passmu, zpos_notmu;
            float trch, tpmu, tpel, tnod, ttid_passmu, trid_passmu;
            double trackdir[3];
            //track labels - reset as 0 before each PID loop
            int mulabel, notmulabel;
            //momentum components for munonCCQE sum
            float xmom, ymom, zmom, psum;
            //(re)Set track type labels = 0
            mulabel    = 0;
            notmulabel = 0;
            //reset momentum components ready to sum over contents of next bunch
            xmom = 0.0; ymom = 0.0; zmom = 0.0; psum = 0.0;
            float Emutracksum = 0.0, Eothertracksum = 0.0;
            
            
            //*#CUT#* VETO
            //Skip the bunch if there is activity in the P0D, P0DECal or TPC1
            if(pod[bunch]   != 0){continue;}
            EventsPassedPODVeto++;
            if(EventID==bob){std::cout << "Marker 1" << endl;}
            
            if(tpc1a[bunch]  != 0){continue;}
            EventsPassedTPC1Veto++;
            if(EventID==bob){std::cout << "Marker 2" << endl;}

            
            
            //PID loop
            for(int j = 0; j < NPIDs; ++j){                            //   ***IN TRACK (PID) LOOP***
                
             //Sets up  a pointer that allows us to access the variables more easily.
                ND::TGlobalReconModule::TGlobalPID *trTrack = (ND::TGlobalReconModule::TGlobalPID*)(*PIDs)[j];
                if(EventID==bob){std::cout << "ID: " << trTrack->UniqueID << " NTPCS: " << trTrack->NTPCs << " Momentum: " << trTrack->FrontMomentum << " bunch: " << buntm.bunchno(trTrack->FrontPosition.T()) << std::endl;}
                
                if(!((fabs(trTrack->FrontPosition.T() - bunch1) < 60.) || (fabs(trTrack->FrontPosition.T() - bunch2) < 60.) || (fabs(trTrack->FrontPosition.T() - bunch3) < 60.) || (fabs(trTrack->FrontPosition.T() - bunch4) < 60.) || (fabs(trTrack->FrontPosition.T() - bunch5) < 60.) || (fabs(trTrack->FrontPosition.T() - bunch6) < 60.) || (fabs(trTrack->FrontPosition.T() - bunch7) < 60.) || (fabs(trTrack->FrontPosition.T() - bunch8) < 60.))){continue;}
                
                //*#CUT#* BUNCH TIMING
                if(bunch != buntm.bunchno(trTrack->FrontPosition.T())){continue;}
                EventsPassedBunchTiming++;
//                if(EventID==bob){
//                    std::cout << "Marker 3 ID: " << trTrack->UniqueID << " NTPCS: " << trTrack->NTPCs << " bunch: " << bunch << " buntm: " << buntm.bunchno(trTrack->FrontPosition.T()) << endl;}
                
                //*#CUT#* TPC INFO (if no TPC info for the track go to the next track (next j) in the entry
                if(trTrack->NTPCs == 0){continue;}
                EventsPassedTPCInfo++;
                if(EventID==bob){std::cout << "Marker 4, No. TPCs: " << trTrack->NTPCs << endl;}
                
                //*#CUT#* VERTEX IN FGD1 FID VOL
                
                
                bool passFGD1= FGD1xmin < trTrack->FrontPosition.X() && trTrack->FrontPosition.X() < FGD1xmax && FGD1ymin < trTrack->FrontPosition.Y() && trTrack->FrontPosition.Y() < FGD1ymax && FGD1zmin < trTrack->FrontPosition.Z() && trTrack->FrontPosition.Z() < FGD1zmax;
                if(!passFGD1){continue;} //if not in FGD1, reject
                EventsPassedVertexCut++;
                if(EventID==bob){
                    std::cout << "Marker 5" << endl;}
                
                //*#CUT#* FORWARD-GOING TRACKS ONLY
                bool isfwd = trTrack->isForward; //This returns 1 if the track forward, 0 if not.
                if(!isfwd){continue;} //CUT forward tracks only
                EventsPassedForwardCut++;
                if(EventID==bob){std::cout << "Marker 6" << endl;}
                
                //Define quantities for cuts
                int tpc1, tpc2, tpc3;
                double charge;
                int nodes = 0;
                float pulle, pullm;
                
                //Get tpc component numbers ready for cuts and get values when in tpc2 for mu sel
                TClonesArray * TPCg = trTrack->TPC;
                tpc1 = 0; tpc2 = 0; tpc3 = 0;//set to 0 ready for cut
                for(int ntpc = 0; ntpc < trTrack->NTPCs; ++ntpc){         //****IN TPC LOOP****
                    ND::TGlobalReconModule::TTPCObject *gtpctrack =(ND::TGlobalReconModule::TTPCObject*)(*TPCg)[ntpc];
                    if(gtpctrack->Detector == 1){++tpc1;}
                    if(gtpctrack->Detector == 3){++tpc3;}
                    if(gtpctrack->Detector == 2){++tpc2; //only interested in 2. Get values ready for next cuts
                        charge = gtpctrack->Charge;
                        pullm  = gtpctrack->PullMuon;
                        pulle  = gtpctrack->PullEle;
                        nodes  = gtpctrack->NNodes;
                    }
                }//end of loop over TPCs                                ****END TPC LOOP****
                
                //*#CUT#* TRACK HAS TPC2 COMPONENT
                if(tpc2 == 0){continue;} //if track has no tpc2 comp, reject
                EventsPassedTPC2Comp++;
                if(EventID==bob){std::cout << "Marker 7" << endl;}
                
                //Muon candidate selection setup
                trch = charge;//do charge, pullm and pulle need defining again?
                tpmu = pullm;
                tpel = pulle;
                tnod = nodes;
                
                bool passchrg = trch < 0;
                bool passmuonpid = tpmu < 2.0 && tpmu > -2.0 && tpel < -2.0;
                bool passtrqual = tnod > 18;
                
                //************************************************
                // Muon candidate selection and energy calculation
                //************************************************
                //*#SELECT#* -ve, pull, 18
                if(passchrg && passmuonpid && passtrqual){
                    ++mulabel;
                    EventsPassedMuonCandCut++;
                    if(EventID==bob){std::cout << "Marker 8, bunch: " << buntm.bunchno(trTrack->FrontPosition.T()) << endl;}
                    
                    // If muon candidate calculate Energy of track based on this
                    Emutracksum += TMath::Sqrt(trTrack->FrontMomentum*trTrack->FrontMomentum + mumass*mumass);
                    
                    // Store momentum and position of highest momentum track
                    if(trTrack->FrontMomentum > pmax_passmu){ 
                        pmax_passmu = trTrack->FrontMomentum;
                        xpos_passmu = trTrack->FrontPosition.X();
                        ypos_passmu = trTrack->FrontPosition.Y();
                        zpos_passmu = trTrack->FrontPosition.Z();
                        trackdir[0] = trTrack->FrontDirection.X();
                        trackdir[1] = trTrack->FrontDirection.Y();
                        trackdir[2] = trTrack->FrontDirection.Z();
                        ttid_passmu = trTrack->TrueParticle.Vertex.ID; //This is for the true vertex
                        trid_passmu = trTrack->TrueParticle.ID; //Trajectory ID for the track
                    }
                }
                
                
                //************************************************
                // Pion candidate selection and energy calculation
                //************************************************
                //If not muon candidate...(assume pion for E calc)
                else{
                    ++notmulabel;
                    
                    // Calculate energy of track assuming it is a pion
                    Eothertracksum += TMath::Sqrt(trTrack->FrontMomentum*trTrack->FrontMomentum + pimass*pimass);
                    
                    // Store position and momentum of track for highest momentum track
                    if(trTrack->FrontMomentum > pmax_notmu){ 
                        pmax_notmu = trTrack->FrontMomentum;
                        xpos_notmu = trTrack->FrontPosition.X();
                        ypos_notmu = trTrack->FrontPosition.Y();
                        zpos_notmu = trTrack->FrontPosition.Z();
                    }
                }
                
                //****************************************************************
                // regardless of track type, add up momentum vectors for all tracks
                //****************************************************************
                psum += trTrack->FrontMomentum; //basic sum to be added to masses
                xmom += trTrack->FrontMomentum*trTrack->FrontDirection.X(); //Assumes FrontDir is normalised vector
                ymom += trTrack->FrontMomentum*trTrack->FrontDirection.Y(); //Assumes FrontDir is normalised vector
                zmom += trTrack->FrontMomentum*trTrack->FrontDirection.Z(); //Assumes FrontDir is normalised vector
                
                
            }//end of PID loop                                          ***END PID LOOP***
            
            //**Process bunches depending on contents**
            
            //***********************************************
            // Select only bunches containing muon candidate
            //***********************************************
            if(mulabel == 0){continue;} //if no mu candidates in bunch, go to next bunch
            EventsPassedMuonCut++;
            std::cout << "Event: " << EventID << std::endl;
            if(EventID==bob){std::cout << "Marker 9" << endl;}
            //std::cout << "Event no: " << EventID << endl;
            
            
            //***************************
            //*#SELECT CCQE MU EVENTS#*
            //***************************
            if(mulabel >= 1 && notmulabel == 0){   //if at least 1 mu track and no other tracks in bunch
                ++NCCQEmu;
                //std::cout << "Event: " << EventID << std::endl;
                //Enu calc using pmax_passmu (so if there was >1 CCQEmu in the bunch, will only use highest mom one)
                double costh = (trackdir[0]*nudir[0])+(trackdir[1]*nudir[1])+(trackdir[2]*nudir[2]);; //|1||1|costh
                double pl = pmax_passmu; //lepton momentum
                double el = TMath::Sqrt(pl*pl + ml*ml); //lepton energy
                
                //nu energy
                double enu = (((mn-V)*el) - (0.5*ml*ml) + (mn*V) - (0.5*V*V) + (0.5*((mp*mp)-(mn*mn)))) / (mn - V - el + (pl*costh));
                
                //Now need to calculate the weight using 2011a histograms.  Need to pass the following to the function
                //TClonesArray for true vertices, Number of true vertices, True ID of highest momentum reconstructed muon, then
                //the 4 histograms for the 2011 weighting
                double weight = fluxw.fweight(Vertices,NVtx,ttid_passmu,numu,numub,nue,nueb);
                //double weight = 1.;
                
                EventsPassedCCQE++;
                
                // *** Fill plots ***
                
                // Energy plots
                hEnu->Fill(enu, weight);
                
                hEest_all->Fill(enu, weight);
                
                //fill quadrant spectra based on highest mom mu track
                if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hEestQ1->Fill(enu, weight);} //fine versions
                if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hEestQ3->Fill(enu, weight);}
                if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hEestQ1_b->Fill(enu, weight);} //binned versions
                if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hEestQ3_b->Fill(enu, weight);}
                
                //find missing events
                int vtxflag = 0;
                //double weight = 1.;
                //get truth beampdg values from rooVertex

//                std::cout << "Event: " << EventID << " trid_passmu: " << trid_passmu << std::endl;
                
                for(int i = 0; i < NVtx; ++i){
                    ND::NRooTrackerVtx *rooVertex = (ND::NRooTrackerVtx*)Vtx->At(i);
                    if(bunch != buntm.vtbunchno(rooVertex->TimeInSpill)){continue;}
                    
                    int vertexid = rooVertex->TruthVertexID;//needed for truth matching
                    int beampdg = rooVertex->NuParentPdg;//get parent type
                    bool passpi = beampdg == 211 || beampdg == -211;
                    bool passK = beampdg == 321 || beampdg == -321 || beampdg == 130;
                    bool passmupar = beampdg == 13 || beampdg == -13;
                    if(ttid_passmu == vertexid){vtxflag = 1;}//to check for missing events
                    //fill quad spectra based on true parent
                    if(vertexid == ttid_passmu && passpi){
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecPiQ1->Fill(enu, weight);}//fine
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecPiQ3->Fill(enu, weight);}
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecPiQ1_b->Fill(enu, weight);}//binned
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecPiQ3_b->Fill(enu, weight);}
                    }
                    if(vertexid == ttid_passmu && passK){
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecKQ1->Fill(enu, weight);}//fine
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecKQ3->Fill(enu, weight);}
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecKQ1_b->Fill(enu, weight);}//binned
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecKQ3_b->Fill(enu, weight);}
                    }
                    if(vertexid == ttid_passmu && passmupar){
                        cout << "Here's the muon parent!!" << endl;
                        cout << "true vertex ID " << ttid_passmu << "\t bunch " << bunch << "\t entry " << ient << endl;
                        cout << "Enu " << enu << "\t weight " << weight << endl;
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecMuQ1->Fill(enu, weight);}
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecMuQ3->Fill(enu, weight);}
                    }
                }
                
                if(vtxflag == 0){//print out if event got lost
                    cout << "LOST EVENT" << endl;
                    cout << "true vertex ID " << ttid_passmu << "\t bunch " << bunch << "\t entry " << ient << endl;
                    cout << "Enu " << enu << "\t weight " << weight << endl;
                }
                
                //Compare to truth values to test E_Est
                //first off need to get the true neutrino energy.
                for(int vt = 0; vt < NTVtx; ++vt){
                    ND::TTruthVerticesModule::TTruthVertex *vtxTrue = (ND::TTruthVerticesModule::TTruthVertex*)(*Vertices)[vt];
                    if(vtxTrue->ID != ttid_passmu){continue;}//go to next true vertex if it doesn't equal ttid
                    //get the PDG code for the track - have to match track true ID with a trajectory ID
                    for(int j = 0; j < NTraj; ++j){
                        ND::TTruthTrajectoriesModule::TTruthTrajectory *trjTrue = (ND::TTruthTrajectoriesModule::TTruthTrajectory*)(*Trajectories)[j];
                        
                        if(trjTrue->PDG == 13 && trjTrue->ID == trid_passmu){
                            hEnuVsTrue->Fill(vtxTrue->NeutrinoMomentum.E(),enu);
                            //non-weighted - comp just asks accuracy of enu calc per event, doesn't care weight of event
                            //old position of parent bits
                        }
                        
                    }         
                }//end of truth comp
                
            }
            
            
            
            
            //*******************************
            //*#SELECT nonCCQE MU EVENTS#*
            //*******************************
            
            //unlikely to have 2 mu events and possible that a non-mu passes the cuts but is actually track crap, so ptot is ok...
            
            //150mm position cut
            double zdiff = zpos_passmu - zpos_notmu;
            bool pass150 = TMath::Abs(zdiff) < 150;
            
            if (mulabel >= 1 && notmulabel > 0 && pass150){//1 mu track AND misc track crap present
                ++NnCCQEmu;
                //std::cout << "Event: " << EventID << std::endl;
                //now put summed x,y,z components of momentum into vector so can plot magnitude
                TVector3 vecsum;
                vecsum.SetXYZ(xmom,ymom,zmom);
                double pbunch = vecsum.Mag(); //vecsum
                double pbunchcomp = pbunch + mulabel*mumass + notmulabel*pimass; //vecsum all tracks + masses
                double psumcomp = psum + mulabel*mumass + notmulabel*pimass; //scal sum all tracks + masses
                
                // If there the event in not CCQE calculate the energy of the tracks in the event.
                double Etracksum = Emutracksum + Eothertracksum; // muon track + pion track
                
                double weight = fluxw.fweight(Vertices,NVtx,ttid_passmu,numu,numub,nue,nueb);
                //double weight = 1.;
                
                EventsPassednonCCQE++;
                
                // *** Fill plots ***
                hpbunch->Fill(pbunch, weight); // ?Momentum of the bunch/Non CCQE Energy
                hEest_all->Fill(Etracksum, weight);//change to Etracksum as most 'proper' estimator
                //std::cout << "Event: " << EventID << " filling hEstall, Etracksum " << Etracksum << std::endl;
                
                //fill quadrant spectra based on position of highest mom mu track
                if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hEestQ1->Fill(Etracksum, weight);}//fine
                if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hEestQ3->Fill(Etracksum, weight);}
                if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hEestQ1_b->Fill(Etracksum, weight);}//binned
                if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hEestQ3_b->Fill(Etracksum, weight);}
                
                int vtxflag = 0;//to check for missing events
                //double weight = 1.;
                //get beampdg values from rooVertex
                for(int i = 0; i < NVtx; ++i){
                    ND::NRooTrackerVtx *rooVertex = (ND::NRooTrackerVtx*)Vtx->At(i);
                    if(bunch != buntm.vtbunchno(rooVertex->TimeInSpill)){continue;}
                    int vertexid = rooVertex->TruthVertexID;//needed for truth matching
                    int beampdg = rooVertex->NuParentPdg;//get parent type
                    bool passpi = beampdg == 211 || beampdg == -211;
                    bool passK = beampdg == 321 || beampdg == -321 || beampdg == 130;
                    bool passmupar = beampdg == 13 || beampdg == -13;
                    if(ttid_passmu == vertexid){vtxflag = 1;}//to check for missing events
                    //fill quad spectra based on true parent
                    if(vertexid == ttid_passmu && passpi){
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecPiQ1->Fill(Etracksum, weight);}//fine
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecPiQ3->Fill(Etracksum, weight);}
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecPiQ1_b->Fill(Etracksum, weight);}//binned
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecPiQ3_b->Fill(Etracksum, weight);}
                    }
                    if(vertexid == ttid_passmu && passK){
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecKQ1->Fill(Etracksum, weight);}//fine
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecKQ3->Fill(Etracksum, weight);}
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecKQ1_b->Fill(Etracksum, weight);}//binned
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecKQ3_b->Fill(Etracksum, weight);}
                    }
                    if(vertexid == ttid_passmu && passmupar){
                        cout << "Here's the muon parent!!" << endl;
                        cout << "true vertex ID " << ttid_passmu << "\t bunch " << bunch << "\t entry " << ient << endl;
                        cout << "Etracksum " << Etracksum << "\t weight " << weight << endl;
                        if(xpos_passmu < FGD1xmid && ypos_passmu > FGD1ymid){hreconspecMuQ1->Fill(Etracksum, weight);}
                        if(xpos_passmu > FGD1xmid && ypos_passmu < FGD1ymid){hreconspecMuQ3->Fill(Etracksum, weight);}
                    }
                }//end beampdg bit
                if(vtxflag == 0){
                    cout << "LOST EVENT" << endl;
                    cout << "true vertex ID " << ttid_passmu << "\t bunch " << bunch << "\t entry " << ient << endl;
                    cout << "Etracksum " << Etracksum << "\t weight " << weight << endl;
                }
                
                //Compare to truth values to test E_Est
                //first off need to get the true neutrino energy.
                for(int vt = 0; vt < NTVtx; ++vt){
                    ND::TTruthVerticesModule::TTruthVertex *vtxTrue = (ND::TTruthVerticesModule::TTruthVertex*)(*Vertices)[vt];
                    if(vtxTrue->ID != ttid_passmu){continue;}//go to next true vertex if it doesn't equal ttid
                    //get the PDG code for the track - have to match track true ID with a trajectory ID
                    for(int j = 0; j < NTraj; ++j){
                        ND::TTruthTrajectoriesModule::TTruthTrajectory *trjTrue = (ND::TTruthTrajectoriesModule::TTruthTrajectory*)(*Trajectories)[j];
                        if(trjTrue->PDG == 13 && trjTrue->ID == trid_passmu){
                            hpbunchVsTrue->Fill(vtxTrue->NeutrinoMomentum.E(),pbunch);//test vec mom sum method
                            hpbunchcompVsTrue->Fill(vtxTrue->NeutrinoMomentum.E(),pbunchcomp);//test scalar mom + mass sum method
                            hpsumcompVsTrue->Fill(vtxTrue->NeutrinoMomentum.E(),psumcomp);//test scalar mom + mass sum method
                            hEtracksumVsTrue->Fill(vtxTrue->NeutrinoMomentum.E(),Etracksum);//test scalar mom + mass sum method
                            //don't weight - just true vs recon per event test
                            //old beampdg position
                        }//end comp with true energy
                    }
                }
                
            }
            
        }//end of bunch loop
        
        
    }//end of entries loop
    
    //****************
    // Prints
    //****************
    
    cout << "Events seen:" << EventsSeen << endl;
    cout << "Events passed POD Veto: " << EventsPassedPODVeto << endl;
    cout << "Events passed TPC1 Veto: " << EventsPassedTPC1Veto << endl;
    cout << "Events passed Bunch timing cut: " << EventsPassedBunchTiming << endl;
    cout << "Events passed TPC info cut: " << EventsPassedTPCInfo << endl;
    cout << "Events passed vertex cut: " << EventsPassedVertexCut << endl;
    cout << "Events passed forward cut: " << EventsPassedForwardCut << endl;
    cout << "Events passed TPC2 Comp cut: " << EventsPassedTPC2Comp << endl;
    cout << "Events passed Muon Cand. Cut: " << EventsPassedMuonCandCut << endl;
    cout << "Events passed Muon Cut: " << EventsPassedMuonCut << endl;
    cout << "Events passed CCQE:" << EventsPassedCCQE << endl;
    cout << "Events passed nonCCQE:" << EventsPassednonCCQE << endl;
    cout << "No. of empty bunches: " << EmptyBunch << endl;
    
    
    //******************
    // Draw Histograms
    //******************
    c1->cd();
    hEnu->Draw();
    c2->cd();
    hpbunch->Draw();
    c3->cd();
    hEest_all->Draw();

    c4->cd(); //fine
    hEestQ3->Draw();
    hEestQ1->Draw("SAME");
    legemuQ1Q3->Draw();
    c5->cd();
    hEnuVsTrue->Draw("colz");
    c6->cd();
    hpbunchVsTrue->Draw("colz");
    c7->cd();
    hpbunchcompVsTrue->Draw("colz");
    c8->cd();
    hpsumcompVsTrue->Draw("colz");
    c9->cd();
    hEtracksumVsTrue->Draw("colz");
    c10->cd();//fine
    hreconspecPiQ3->Draw();
    hreconspecPiQ1->Draw("SAME");
    hreconspecKQ3->Draw("SAME");
    hreconspecKQ1->Draw("SAME");
    legePiKQ1Q3->Draw("SAME");
    c11->cd(); //binned
    hEestQ3_b->Draw();
    hEestQ1_b->Draw("SAME");
    c12->cd();//binned
    hreconspecPiQ3_b->Draw();
    hreconspecPiQ1_b->Draw("SAME");
    hreconspecKQ3_b->Draw("SAME");
    hreconspecKQ1_b->Draw("SAME");
    legePiKQ1Q3_b->Draw("SAME");
    
    
    //*******************
    // Save root file
    //*******************
    TFile *file = new TFile("/home/dewhurst/WorkArea/RebekahsAnalysis/OutputFiles/mymacro12out_test.root","RECREATE");
    hEnu->Write();
    hpbunch->Write();
    hEest_all->Write();
    hEestQ1->Write();
    hEestQ3->Write();
    hEestQ1_b->Write();
    hEestQ3_b->Write();
    hEnuVsTrue->Write();
    hpbunchVsTrue->Write();
    hpbunchcompVsTrue->Write();
    hpsumcompVsTrue->Write();
    hEtracksumVsTrue->Write();
    hreconspecPiQ3->Write();
    hreconspecPiQ1->Write();
    hreconspecKQ3->Write();
    hreconspecKQ1->Write();
    hreconspecMuQ3->Write();
    hreconspecMuQ1->Write();
    //binned versions
    hEestQ1_b->Write();
    hEestQ3_b->Write();
    hreconspecPiQ3_b->Write();
    hreconspecPiQ1_b->Write();
    hreconspecKQ3_b->Write();
    hreconspecKQ1_b->Write();
    
    file->Close();
    

    

    
    


    
    

    
}

//Done!
